Per dataset che hanno un numero piccolo o moderato di predittori, è possibile valutare tutte le coppie di interazioni.
Tuttavia maggiore è il numero di termini di interazione che vengono valutati, maggiore è la probabilità di trovare delle interazioni dette false positive, ovvero che risultano statisticamente significative ma non a causa di una vera relazione e possono dunque peggiorare le performance predittive del modello.

Simple Screening

Gli approcci tradizionali per individuare i termini di interazioni più rilevanti si basano sui modelli nidificati.
Si consideri un modello lineare con 2 soli regressori, si definisce main effects model il seguente modello:
\(M_1: y=\beta_0+\beta_1x_1+\beta_2x_2+\epsilon\)

E’ possibile costruire un secondo modello che aggiunge una potenziale interazione tra i due regressori.
\(M_2: y=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_1x_2+\epsilon\)

Questi due modelli sono detti “nidificati” perché il primo modello è incluso nel secondo (\(M_1 \subseteq M_2\)).
In presenza di questa struttura, è possibile effettuare un confronto statistico sulla quantità di informazione addizionale che è stata acquisita dal termine di interazione.
In questo contesto, con riferimento a una regressione lineare, è necessario confrontare l’errore residuo di \(M_1\) e \(M_2\) per valutare se il miglioramento dell’errore,ponderato per i gradi di libertà, è sufficiente per essere considerato reale.
Per la regressione lineare, la funzione obiettivo utilizzata per confrontare i modelli è la verosimiglianza statistica (in questo caso, l’errore residuo).
Per altri modelli, ad esempio il modello di regressione logistica, la funzione obiettivo per confrontare i modelli nidificati coincide con la verosimiglianza binomiale.
Il risultato di questo test statistico lo si legge dal p-value che può essere interpretato come il tasso di risultati falsi positivi e riflette la probabilità che l’informazione aggiuntiva catturata dal termine di interazione sia dovuta alla casualità.
Più è basso il p-value minore è la probabilità che l’informazione aggiuntiva catturata dal termine di interazione sia dovuta alla casualità.
Una evoluzione del modello nidificato sfrutta la cross-validation in modo da creare più coppie di modelli annidati a partire da differenti versioni del training set, a differenza dell’approccio tradizionale in cui gli stessi dati usati per la creazione del modello vengono riutilizzati per la sua valutazione che può essere svolta mediante qualsiasi misura di performance.
Considerato che la valutazione dipende molto dal contesto, questi due approcci non garantiscono conformità nei risultati.
Come detto in precedenza più sono i confronti che vengono effettuati, più alta è la possibilità di trovare interazioni false-positive.
Esistono diversi metodi per affrontare questo problema, ad un estremo, si sceglie di non effettuare controlli per i falsi positivi, oppure si applica la correzione di Bonferroni che utilizza una “severa” penalità esponenziale per minimizzare eventuali risultati falsi positivi.
Tuttavia la tecnica più equilibrata per individuare e gestire i falsi positivi si basa sul false discovery rate (FDR).

Esempio dataset Ames

Nella presente applicazione vengono considerati solamente 18 regressori presenti nel dataset Ames, rispettivamente: Bldg_Type, Neighborhood, Year_Built, Gr_Liv_Area, Full_Bath, Year_Sold, Lot_Area, Central_Air, Longitude, Latitude, MS_SubClass, Alley, Lot_Frontage, Pool_Area, Garage_Finish, Foundation, Land_Contour e Roof_Style.
Il comando recipe() permette di applicare opportune trasformazioni ai dati iniziali.

ames_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05)  %>%
  step_bs(Longitude, Latitude, options = list(df = 5)) %>% 
  step_zv(all_predictors())

In particolare, le trasformazioni effettuate sono state le seguenti:

## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         18
## 
## Operations:
## 
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood
## B-Splines on Longitude, Latitude
## Zero variance filter on all_predictors()

Si procede creando una dataframe contenente tutte le possibile coppie di interazioni \(\frac{p(p-1)}{2}\) tra i regressori.

ames_preds <- ames_rec$var_info$variable[ames_rec$var_info$role == "predictor"]
interactions <- t(combn(ames_preds, 2))
colnames(interactions) <- c("var1", "var2")

Si decide di addestrare un modello lineare sui dati del training set utilizzando una repeated k-fold cross-validation, iterando la procedura 5 volte con K=10.

int_ctrl <- trainControl(method = "repeatedcv", repeats = 5)
set.seed(2691)
main_eff <- train(ames_rec, 
                  data = ames_train, 
                  method = "lm", 
                  metric = "RMSE",
                  trControl = int_ctrl)

Si implementa la funzione compare_models_1way al fine di confrontare il modello senza interazioni (main effect model) e il main effect model più un’ interazione.

compare_models_1way <- function(a, b, metric = a$metric[1], ...) {
  mods <- list(a, b)
  rs <- resamples(mods)
  diffs <- diff(rs, metric = metric[1], ...)
  diffs$statistics[[1]][[1]]
}

In questo modo è possibile valutare l’effetto delle interazioni in termini di: Reduction, Pvalue, RMSE, anova_p.

for (i in 1:nrow(interactions)) {
  tmp_vars <- c("Class", interactions$var1[i], interactions$var2[i])
  
  int <-
    as.formula(paste0(
      "~ starts_with('",
      interactions$var1[i],
      "'):starts_with('",
      interactions$var2[i],
      "')"
    ))
  
  int_rec <- 
    ames_rec %>% 
    step_interact(int) %>% 
    step_zv(all_predictors())
  
  set.seed(2691)
  main_int <- train(int_rec, 
                    data = ames_train, 
                    method = "lm", 
                    metric = "RMSE",
                    trControl = int_ctrl)  
  
  tmp_diff <- compare_models_1way(main_int, main_eff, alternative = "less")
  interactions$RMSE[i] <- getTrainPerf(main_eff)[1, "TrainRMSE"]
  interactions$Reduction[i] <- -tmp_diff$estimate
  interactions$Pvalue[i] <- tmp_diff$p.value
  
  a1 <- 
    anova(main_eff$finalModel, main_int$finalModel) %>% 
    tidy() %>% 
    slice(2) %>% 
    pull(p.value)
  
  interactions$anova_p[i] <- a1
}

head(interactions)

## # A tibble: 6 × 6
##   var1      var2          Reduction      Pvalue   RMSE  anova_p
##   <chr>     <chr>             <dbl>       <dbl>  <dbl>    <dbl>
## 1 Bldg_Type Neighborhood  0.000184  0.0479      0.0748 2.08e- 4
## 2 Bldg_Type Year_Built   -0.0000263 0.680       0.0748 7.94e- 2
## 3 Bldg_Type Gr_Liv_Area   0.000736  0.000000471 0.0748 7.14e-10
## 4 Bldg_Type Full_Bath     0.000215  0.00107     0.0748 1.53e- 3
## 5 Bldg_Type Year_Sold    -0.0000560 0.915       0.0748 2.33e- 1
## 6 Bldg_Type Lot_Area     -0.0000249 0.683       0.0748 6.88e- 2

plot <- 
  ggplot(all_p, aes(x = Value)) + 
  geom_histogram(breaks = (0:20)/20) + 
  geom_rug(col = "red") + 
  facet_grid(Estimate ~ Method)

Questa rappresentazione confronta la distrubuzione dei p-value ottenuta con l’approccio tradizionale e la cross-validation, valutando le diverse tipologie di aggiustamento (No Adjustment, FDR, Bonferroni).
Emerge che utilizzando l’approccio tradizionale quasi tutti gli effetti delle interazioni vengono valutate significative, ciò potrebbe essere dovuto a una potenziale presenza di overfitting. L’approccio di cross validation ha un tasso minore di significatività ed è meno influenzato da entrambi i tipi di aggiustamenti.

Valutare i termini di interazione singolarmente impedisce ai modelli semplici di valutare l’importanza dei termini di interazione nel modello completo ( contenente tutti i main effects e i relativi termini di interazione), per questo se la completa enumerazione è possibile si potrebbe provare un approccio che ricerca le interazioni partendo da un modello completo.

Utilizzando quest’approccio potrebbe però capitare che nei dati si presentino più predittori rispetto alle osservazioni del campione, ciò potrebbe non essere un problema per modelli come basati su alberi decisionali, neural networks, SVM, e K-nearest neighbors, che sono applicabili anche quando i dati contengono più predittori rispetto alla numerosità delle osservazioni del campione. Tuttavia altre tecniche più interpretabili e che spesso riportano buone performance, ad esempio le regressioni lineari e logistiche, non possono essere utilizzate direttamente sotto queste condizioni. Questa limitazione deriva dal fatto che qualora si hanno più predittori che osservazioni nel campione, oppure un predittore può essere scritto come combinazione di altri, non è possibile effettuare l’inversione della matrice del disegno.

Essendo l’interpretabilità e la bontà del modello obiettivi della feature engineering, si sono cercati metodi per utilizzare anche in questo caso la regressione lineare e logistica. In questo contesto, si introduce una famiglia dei tecniche di modellazione dette modelli penalizzati.

I metodi più utilizzati tra i modelli penalizzati sono la ridge regression e la lasso regression. Questi due tipi differenti di penalità sono solitamente associati a differenti task:

Ridge penality, utilizzata per combattere la collinearità tra predittori;
Lasso penality, utilizzata per l’eliminazione dei predittori superflui.

Talvolta potrebbe essere necessario combinare le due penalità, a tal fine è stato introdotto il modello glmnet.

Glmnet regression

Nel modello glmnet compaiono entrambi i termini di penalizzazione ridge e lasso, inoltre, sono presenti due parametri di tuning:

\(\lambda\) la penalità totale (\(\lambda\)=\(\lambda_r\)+\(\lambda_l\))
\(\alpha\) che indica la proporzione di lambda che viene associata alla regressione lasso. \((1-\alpha)\) è la proporzione di lambda associata alla regressione ridge.

Ad esempio scegliendo \(\alpha=1\) avremo un modello fully lasso penality, scegliendo \(\alpha=0.5\) avremo un modello che è metà lasso penality e metà ridge penality.

La regressione tramite modello glmnet si può vedere, come nel caso della regressione lineare, come un problema di minimizzazione, in particolare del valore del somma dei quadrati dei residui, calcolata con la seguente equazione:

\[ SSE=\sum_{i=1}^n (y_i - \hat{y}_i)^2 +(1-\alpha)\lambda_r \sum_{j=1}^P \beta^2_j + \alpha \lambda_l\sum_{j=1}^P |\beta_j| \]

Esempio con dati ames

L’esempio è stato svolto utilizzando i dati ames e consiste nel tuning di due modelli glmnet. Un primo modello contenente solo i main effect e un secondo modello contenente main effect e relative interazioni. Il modello glmnet dopo aver provato tutte le combinazioni di valori dei parametri di tuning \(\alpha\) e \(\lambda\) sceglie il “best tuning”, cioè il modello che minmizza il valore dell’RMSE. La scelta dei parametri di tuning ottimali permette di selezionare le variabili esplicative più significative e di stimare i loro coefficienti.

Recupero dati e creazione modello main

Il primo passo è quello di recuperare i dati ames e di effettuare lo split in training set e test set e successivamente con il comando vfold_cv() in 10 folds per poi svolgere una cross validation.

Si passa successivamente alla creazione del primo modello, come detto in precedenza, creato utilizzando solo i main effect tramite una struttura chiamata recipe() che permette inoltre di specificare una pipeline per il preprocessing dei dati. In questo passo si specificano quindi la variabile target e quelle esplicative del modello e le operazioni da effettuare in fase di preprocessing.

main_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>%
  step_zv(all_predictors()) %>%
  step_bs(Longitude, Latitude, options = list(df = 5)) %>%
  step_center(all_predictors()) %>%
  step_scale(all_predictors())

Creazione interazioni e modelli main + interazioni

Dopo aver creato un dataframe con tutte le possibili interazioni di secondo ordine tra i main effects del primo modello, si crea una nuova recipe() per il secondo modello.

int_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>%
  step_interact(interactions) %>% # Aggiunta delle interazioni
  step_zv(all_predictors()) %>%
  step_bs(Longitude, Latitude, options = list(df = 5)) %>%
  step_center(all_predictors()) %>%
  step_scale(all_predictors())

Tuning dei modelli

Il modello glmnet, utilizza due diversi parametri di tuning \(\alpha\) e \(\lambda\), si crea un griglia contenente tutte le possibili combinazioni di questi parametri che il modello utilizzerà per effettuare il tuning ricercando la coppia che minimizza il valore dell’RMSE. Il parametro \(\alpha\) può assumere valori da \(0.2\) a \(1\), con passo di \(0.2\). Mentre \(\lambda\) valori da \(10^{-4}\) a \(10^{-1}\) con il valore dell’esponente che varia di \(0.1\).

glmn_grid <- expand.grid(alpha = seq(.2, 1, by = .2), lambda = 10^seq(-4, -1, by = 0.1))

Modello main

main_glmn <- 
  train(main_rec, # La recipe 
        data = ames_train, # Dati da utilizzare per il tuning dei parametri
        method = "glmnet", 
        tuneGrid = glmn_grid, # Dataframe con le combinazioni dei parametri di tuning
        trControl = ctrl # Divisione in fold per la cross validation
  )

## $all
## [1] 58
## 
## $main
## [1] 42
## 
## $perf
##   TrainRMSE TrainRsquared   TrainMAE method
## 1 0.0754791     0.8184643 0.05187811 glmnet

Dopo aver effettuato il tuning dei parametri per il primo modello, il modello risultante che minimizza il valore dell’RMSE è composto da 42 main effects che sono stati selezionati tramite la penalizzazione lasso.

Modello main + interazioni

int_glmn <- 
  train(int_rec,
        data = ames_train, 
        method = "glmnet",
        tuneGrid = glmn_grid,
        trControl = ctrl
  )

## $all
## [1] 981
## 
## $main
## [1] 9
## 
## $int
## [1] 132
## 
## $perf
##    TrainRMSE TrainRsquared   TrainMAE method
## 1 0.07589794      0.817775 0.05157754 glmnet

Il secondo modello ottimale selezionato tramite i parametri di tuning è composto invece da 9 main effects e 132 termini di interazione.

Risultato main + interazioni

Vengono mostrati i parametri che permettono di minimizzare il valore dell’RSME.

int_glmn$bestTune

##    alpha      lambda
## 50   0.4 0.006309573

Il plot che mostra il rapporto tra il valore dell’RMSE (asse y) e quello di \(\alpha\) (colori) e \(\lambda\) (asse x).

Creazione modello two-way

Dopo aver effettuato il tuning dei precedenti due modelli e trovato per ogni modello i regressori significativi, l’esempio procede con un modello glmnet two way. Per la sua creazione vengono presi i main effects risultati significativi nel primo modello, si creano tutte le possibili interazioni binarie e infine si aggiungono ad un modello contenente i main effect. Si riesegue poi il tuning su questo terzo modello.

two_stage_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>%
  step_interact(interaction_subset) %>% 
  step_zv(all_predictors()) %>%
  step_bs(Longitude, Latitude, options = list(df = 5)) %>%
  step_center(all_predictors()) %>%
  step_scale(all_predictors())

## $all
## [1] 741
## 
## $main
## [1] 9
## 
## $int
## [1] 108
## 
## $perf
##    TrainRMSE TrainRsquared   TrainMAE method
## 1 0.07614744     0.8166224 0.05192411 glmnet

Il numero di interazioni da studiare cresce esponenzialmente al crescere del numero di variabili esplicative rendendo l’esplorazione dei termini di interazione molto onerosa dal punto di vista computazionale.

Dal lato statistico non sarebbe una scelta saggia studiare tutte le possibili interazioni. Le interazioni importanti tra le variabili esplicative, infatti, sono generalmente poche e all’aumentare del numero di interazioni irrilevanti esaminate aumenta la penalità imposta al fine di valutarle. Di conseguenza, interazioni realmente significative potrebbero essere perse. In questi casi sono necessari altri approcci per rendere la ricerca più efficiente, ma comunque efficace, nel trovare le interazioni importanti.

I metodi presentati finora sono efficaci per individuare interazioni in presenza di un numero moderato di regressori. Metodi basati sugli alberi decisionali risultano particolarmente efficaci per scoprire potenziali interazioni in presenza di molti regressori.
Queste tecniche si basano sul partizionamento ricorsivo che consiste nell’individuare il punto di split ottimo di un predittore che separa le osservazioni in gruppi omogenei rispetto alla variabile risposta.
Una volta effettuato lo split, la stessa procedura si ripete per ciascun sottoinsieme di osservazioni in corrispondenza di ogni nodo successivo.
Questo processo termina quando si raggiunge il numero minimo di osservazioni in ciascun nodo oppure viene soddisfatto un criterio di ottimizzazione.
Quando i dati vengono splittati ricorsivamente, ogni nodo successivo può essere interpretato come una interazione locale tra il nodo precedente e il nodo corrente.
Inoltre, più è frequente il verificarsi di nodi successivi della stessa coppia di predittori, maggiore è la probabilità che l’interazione tra questi sia di tipo globale.
Formalmente, al fine di comprendere l’efficacia della tecnica presentata per individuare una relazione sinergica tra due regressori, si consideri un esempio in cui l’albero di partizionamento ricorsivo ottimo (Figura 1) risulta:

Figura 1: Albero di partizione ricorsiva utile ad individuare una relazione sinergetica tra due predittori

L’interazione viene catturata attraverso diversi splits alternati tra i due predittori ai livelli successivi dell’albero.
Considerando il loro tipo di equazioni, i modelli basati sugli alberi decisionali possono essere interpretati come una interazione pura.
Ad esempio, l’equazione del nodo 4 può essere scritta come:
\[node_4= I(x_1 <0.655)*I(x_2<0.6)*I(x_0.316)*0.939 \]
dove 0.939 coincide con la media delle osservazioni presenti nel nodo 4 mentre \(I\) rappresenta la funzione indicatrice.

Una relazione sinergica produce delle curve di livello non lineari per la variabile risposta (Figura 2), il modello di partizionamento ricorsivo rompe lo spazio in regioni rettangolari, pertanto risultano necessarie più regioni al fine di spiegare correttamente la risposta, questo però limita fortemente la capacità di modellare interazioni globali.

Figura 2 : Risposta prevista ottenibile a partire da un modello di partizione ricorsiva tra 2 regressori che sono in relazione sinergica

Questo grafico rappresenta le previsioni di ciascun nodo terminale rappresentate da regioni rettangolari in cui il colore indica il valore della risposta prevista per la corrispondente regione a cui sono sovrapposte le curve di livello generate dal modello lineare utilizzato per stimare il termine di interazione.
Nonostante il modello stimato, decomponendo lo spazio della variabili in regioni rettangolari, approssimi adeguatamente la relazione sinergica, sono richieste troppe regioni per spiegare questa relazione.
Un modello basato su un singolo albero decisionale può essere caratterizzato da un’alta variabilità, ciò rende le previsioni poco affidabili.
Per superare questo limite, è possibile ricorrere a metodi di ensemble learning come bagging e boosting che riescono ad identificare le principali relazioni tra i predittori, garantendo una migliore capacità previsiva della variabile risposta.
Le tecniche di ensemble learning sono particolarmente efficaci nell’identificare le relazioni tra predittori e risposta perché aggregano molti alberi decisionali utilizzando delle versioni leggermente diverse dei dati originali.
Considerata una relazione sinergica tra x1 e x2, la Figura 3 mostra che i modelli boosted tree e bagged tree sono in grado di approssimare meglio la relazione tra i predittori e la variabile risposta rispetto a un singolo albero decisionale.

Figura 3 : Modello boosted tree vs. Modello bagged tree nel caso di una relazione sinergica tra due predittori

Infatti, la radice quadrata dell’errore quadratico medio (RMSE) per questi 3 modelli risulta:
\(RMSE_{boosting}=0.27\)
\(RMSE_{bagging}=0.6\)
\(RMSE_{tree}=0.9\)

Nonostante la loro complessità e le loro difficoltà nell’identificare le interazioni globali, quando si verificano interazioni importanti all’interno di uno o più sottoinsiemi dei campioni i metodi basati sugli alberi decisionali risultano particolarmente adeguati per modellare questo tipo di interazioni (interazioni locali).
In conclusione, gli alberi decisionali e i metodi di ensemble learning permettono di ottenere delle informazioni rilevanti riguardo le principali interazioni tra regressori, è suggeribile utilizzare tali informazioni in modo da costruire un modello lineare più semplice.
Un altro approccio basato sui metodi di ensemble learning introdotto da Friedman & Popescu (2008) utilizza la dipendenza parziale al fine di identificare le principali interazioni.
In altre parole, questa tecnica confronta l’effetto congiunto di due (o più) predittori con l’effetto individuale di ciascun predittore di un modello.
Se un predittore non interagisce con nessun altro allora la differenza tra l’effetto congiunto e quello individuale sarà vicina a 0.
Qualora, invece, fosse presente un interazione tra un predittore ed un altro allora questa differenza sarebbe maggiore di 0.
Questo confronto avviene per mezzo della statistica H.
Questa statistica presenta distribuzione bimodale qualora si è in presenza di valori sparsi, e solitamente, ciò accade quando un predittore non è coinvolto in nessuna interazione.
Quando il dataset a disposizione è relativamente piccolo, si calcola la statistica H come la mediana delle statistiche H relative a più campioni di tipo bootstrap.
Successivamente viene fissato un cut off relativamente alla statistica H in modo da selezionare le variabili.
Le variabili a cui corrisponde un valore della statistica H maggiore del cut off vengono segnalate come coinvolte in una potenziale interazione.
In questo modo, è possibile inserire in un modello lineare tutte le coppie di interazioni segnalate.

Esempio dati ames

Per comprendere la metodologia presentata si consideri la seguente applicazione:

h_stat <- function(data, ...) {
  require(foreach)
  pre_obj <- pre(data = data, ...)
  
  pred_names <- pre_obj$x_names
  pred_df <- 
    tibble(
      Predictor = pred_names,
      H = 0
    )
  
  int_obj <-
    try(
      interact(
        pre_obj,
        parallel = TRUE,
        plot = FALSE
      ),
      silent = TRUE)
  
  if (!inherits(int_obj, "try-error")) {
    res <-
      tibble(
        Predictor = names(int_obj),
        H = as.numeric(int_obj)
      )
    missing_pred <- setdiff(pred_names, res$Predictor)
    if (length(missing_pred) > 0) {
      res <-
        pred_df %>% 
        dplyr::filter(Predictor %in% missing_pred) %>%
        bind_rows(res)
    }
  } else {
    print(int_obj)
    res <- pred_df
  }
  res
}

boot_h <- function(split, ...) {
  dat <-
    recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
             Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
             Central_Air + Longitude + Latitude + MS_SubClass +
             Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
             Foundation + Land_Contour + Roof_Style,
           data = analysis(split)) %>%
    step_log(Sale_Price, base = 10) %>%
    step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
    step_other(Neighborhood, MS_SubClass, Roof_Style, Foundation,
               threshold = 0.05)  %>%
    step_zv(all_predictors()) %>% 
    prep(training = analysis(split), retain = TRUE) %>%
    juice()
  h_stat(data = dat, ...)
}

Per motivi computazionali, in questa applicazione vengono considerati solamente 18 regressori presenti nel dataset Ames, rispettivamente: Bldg_Type, Neighborhood, Year_Built, Gr_Liv_Area, Full_Bath, Year_Sold, Lot_Area, Central_Air, Longitude, Latitude, MS_SubClass, Alley, Lot_Frontage, Pool_Area, Garage_Finish, Foundation, Land_Contour e Roof_Style.

ames <- make_ames()
set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)

ames_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, MS_SubClass, Roof_Style, Foundation, 
             threshold = 0.05)  %>%
  step_zv(all_predictors())

ames_pre <- ames_rec %>% prep(ames_train) %>% juice()
ames_rec

## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         18
## 
## Operations:
## 
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood, MS_SubClass, Roof_Style, Foundation
## Zero variance filter on all_predictors()

Il comando recipe() permette di applicare opportune trasformazioni ai dati iniziali.

Dopo aver trasformato le variabili si calcola la statistica H sul training set.

set.seed(3108)  
 obs_h <-  
   h_stat(  
     data = ames_pre,  
     Sale_Price ~ .,  
     normalize = TRUE,  
     ntrees = 10,  
     winsfrac = 0,  
     verbose = TRUE,  
     par.final = TRUE  
   )

Considerato che il dataset a disposizione è relativamente piccolo, si procede calcolando la statistica H anche come la mediana delle statistiche H relative a più campioni di tipo bootstrap.

h_stats <- 
  h_vals %>%
  pull(stats) %>%
  bind_rows() %>%
  group_by(Predictor) %>%
  summarise(Bootstrap = median(H)) %>%
  full_join(
    obs_h %>% dplyr::rename(`Original Data` = H)
  ) %>%
  arrange(Bootstrap)

ordered_pred <- 
  h_stats %>%
  pull(Predictor)

h_stats <-
  h_stats %>%
  gather(Estimator, H, -Predictor) %>%
  mutate(Predictor = factor(Predictor, levels = ordered_pred)) %>% 
  dplyr::filter(H > 0)

Successivamente, viene fissato un cut off pari a 0.001 in modo da segnalare le variabili che risultano coinvolte in una potenziale interazione.

p <- 
  h_stats %>% 
  mutate(
    Predictor = clean_value(Predictor),
    Predictor = str_replace(Predictor, ":", "")
  ) %>% 
  ggplot(aes(x = reorder(Predictor, H), y = H, shape = Estimator)) + 
  geom_hline(yintercept = 0.001, col = "red", alpha = 0.5, lty = 2) + 
  geom_point() + 
  coord_flip() + 
  xlab("") + 
  scale_shape_manual(values = c(16, 1)) + 
  theme(legend.position = "top")
p

Da questa rappresentazione grafica emerge che 6 predittori sono coinvolti in almeno una interazione.

Un secondo approccio utile ad identificare le interazioni combina le informazioni riguardo l’importanza dei predittori con i principi di gerarchia, sparsità ed ereditarietà. Se una interazione viene considerata importante, allora un metodo basato sugli alberi di decisione utilizzerà i predittori coinvolti più volte in più alberi in modo da scoprire la loro relazione con la variabile risposta.
I coefficienti di regressione associati a tali predittori risulteranno significativamente diversi da 0.
In questo modo è possibile creare le principali coppie di interazioni e valutare la loro rilevanza nelle previsione della variabile risposta.
Ogni metodo basato sugli alberi decisionali è caratterizzato da uno specifico algoritmo utile a calcolare l’importanza di un predittore.
In questo contesto, uno degli algoritmi più popolari si ottiene a partire da un random forest che utilizza campioni bootstrap per combinare le previsioni ottenute da diversi modelli, riuscendo così a migliorare la capacità previsiva del modello.
Un albero decisionale costruito su un campione bootstrap utilizza \(\frac{2}{3}\) delle osservazioni originali.
E’ possibile prevedere la risposta della i-esima osservazione utilizzando gli alberi per i quali questa risulta Out-of-Bag.
Questo produce circa \(\frac{B}{3}\) previsioni per la i-esima osservazione che vengono aggregate calcolando la media oppure la moda.
Per valutare l’importanza di uno specifico predittore, i suoi valori vengono mescolati, e successivamente viene ricalcolato l’errore Out-of-Bag.
Qualora il predittore sia rilevante, tale errore dovrebbe peggiorare in modo significativo.
La random forest conduce questa analisi per ogni predittore presente nel modello e crea uno score per ciascuna feature, all’aumentare dello score aumenta il grado di importanza del predittore.

Innanzitutto, si procede stimando un modello di tipo random forest e si ottiene uno score che riflette l’importanza di ciascun regressore.

rf_mod <- ranger(Sale_Price ~ ., data = ames_pre, num.trees = 2000,
                 importance = "impurity", num.threads = workers,
                 seed = 8516)

rf_imp <- 
  tibble(
    Predictor = names(rf_mod$variable.importance),
    Importance = unname(rf_mod$variable.importance)
  )

p <- 
  ggplot(rf_imp, aes(x = reorder(Predictor, Importance), y = Importance)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  xlab("")
p

Questa figura mostra i 10 regressori che risultano più importanti, alcuni di loro sono stati segnalati come potenziali termini di interazione da altre metodologie presentate.
I predittori più importanti sono living area e year built, è ragionevole assumure che una possibile interazione tra questi regressori risulti molto significativa.
In generale, le interazioni tra i principali predittori potrebbero essere incluse in un modello più semplice in modo da valutare se migliora la capacità previsiva della variabile risposta.

Il problema della selezione delle interazioni più importanti, anche considerato un numero moderato di features, risulta impossibile da affrontare con una valutazione completa di tutti i possibili modelli.

Quando si utilizzano modelli lineari si sceglie spesso di utilizzare metodi di forward, backward e stepwise selection.

Miller ha proposto una tecnica alternativa per risolvere i problemi del modello lineare. Anziché rimuovere un predittore, questa tecnica utilizza un approccio sostitutivo: ogni predittore selezionato è scelto sistematicamente per sostituirne un altro all’interno del modello.

Al fine di comprendere la tecnica presenta viene riportato un esempio in cui partendo da un modello formato da dieci variabili esplicative si cerca il miglior modello con solo tre predittori. Un’iterazione dell’algoritmo si può schematizzare come segue:

Si parte selezionando tre variabili esplicative casualmente dalle dieci iniziali e si calcola una misura di performance relativa al modello costruito utilizzandole.
Dopodiché si fissano le prime due esplicative e si sostituisce la terza a rotazione con ciascuna delle rimanenti variabili, calcolando la misura di performance relativa ad ogni modello.
Se nessuno di questi modelli fa meglio del primo, allora la terza variabile esplicativa è quella della prima iterazione, altrimenti è la variabile contenuta nel modello con migliori performance.
Nello step successivo si scarta una delle due variabili precedentemente bloccate (ad esempio la prima) e si bloccano le altre due ripetendo la procedura precedente.
Questo processo viene applicato fintantoché l’ultima variabile del modello è stata selezionata.

Non è garantito che l’ottimo trovato sia un ottimo globale, pertanto bisognerebbe ripetere la procedura con diverse partenze casuali, cioè con combinazioni di variabili iniziali diverse e determinare la soluzione ottima raggiunta.

Hawkins ha esteso questo algoritmo, creando un approccio di modellazione chiamato Feasible Solution Algorithm (FSA), uno dei principali vantaggi è la dimensione dello spazio di ricerca nell’ordine \(q*m*p\) , molto più piccolo dell’intero spazio \(p^m\).

q = numero di partenze casuali
m = numero di termini nel subset
p = numero di variabili esplicative

L’algoritmo precedentemente presentato non considera le interazioni, Lambert lo ha generalizzato ampliando lo spazio di ricerca anche ad esse. L’approccio inizia con un modello base che include le variabili ritenute significative per la previsione della variabile target e, come prima cosa, viene scelto l’ordine delle interazioni tra i main effects. Il processo di identificazione delle interazioni segue la logica dell’FSA orginale, viene riportato il seguente esempio:

Dato un modello base che include tre diverse variabili esplicative si pone l’obiettivo di identificare le interazioni tra coppie di regressori.
L’algoritmo di Lambert seleziona casualmente due regressori e calcola la performance del modello che include i termini di base e il termine di interazione.
Fissato il primo predittore che compone l’interazione, si sostituisce sistematicamente il secondo con i restanti calcolando le corrispondenti performance.
Se il modello originale risulta il migliore si mantiene il primo termine di interazione, altrimenti lo si sostituisce con quello associato al modello migliore.
Si ripete poi il procedimento sostituendo sistematicamente un predittore diverso, fino a convergere in un ottimo.
Si utilizza infine una differente partenza casuale e si riapplica il processo.

Le soluzioni/gli ottimi trovati con le diverse partenze casuali vengono utilizzati per identificare le potenziali soluzioni.

Esempio con dati ames

Anche per l’algoritmo FSA l’esempio viene effettuato sfruttando i dati ames. Per poter confrontare i risultati ottenuti con quelli delle procedure precedenti si è scelto di trasformare anche in questo caso le variabili nominali in variabili dummy individuali e poi procedere alla creazione delle interazioni (la scelta delle variabili dummy, per la creazione delle interazioni, è stata vincolata in modo che non possano essere accoppiate variabili derivanti dallo stesso predittore). Nell’implementazione dell’algoritmo sono state utilizzate 50 partenze casuali, con un subset di massimo 30 variabili esplicative, questo porta ad avere in ogni partenza casuale un massimo di 60 swaps delle variabili che compongono l’interazione.

Recupero dati e creazione modello base

Il primo passo consiste nel recuperare i dati ames ed effettuare lo split in training set e test set e successivamente con il comando vfold_cv() in 10 folds per poi svolgere una cross validation.

Tramite la struttura recipe() si crea la matrice delle variabili esplicative relativa al modello contenente solo main effect. Grazie allo stesso comando si crea la pipeline di preprocessing dei dati, in particolare le variabili nominali sono trasformate in variabili dummies individuali.

ames_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built + 
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>% 
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>% 
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>% 
  step_bs(Longitude, Latitude, options = list(df = 5))  %>% 
  step_zv(all_predictors())

ames_rec

## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         18
## 
## Operations:
## 
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood
## Dummy variables from all_nominal()
## B-Splines on Longitude, Latitude
## Zero variance filter on all_predictors()

Viene mostrato il contenuto della recipe(), in questo caso il modello è formato da una variabile target e 18 predittori.

Inizializzazione parametri

Con il seguente comando viene creato un set di metriche che verranno poi utilizzate nell’esecuzione dell’algoritmo.

multi_metric <- metric_set(rmse, mae, rsq)

Un’altra informaznione necessaria per l’algoritmo è la specificazione del modello di regressione da utilizzare. Nel nostro caso è stata usata una regressione lineare.

lr_spec <- linear_reg() %>% 
  set_engine("lm")

Esecuzione dell’algoritmo fsa

Si esegue l’algoritmo fsa appoggiandosi alla funzione fsa_two_way contenuta nel file fsa_function. La funzione esegue un fsa two way, in quanto si cercano solo le interazioni di secondo grado tra i predittori. Per l’esecuzione della suddetta funzione vengono forniti i seguenti parametri:

ames_folds = la divisione in fold per la cross validation
ames_rec = la pipeline per il preprocessing
lr_spec = il modello da utilizzare per la linear regression
multi_metric = le metriche di valutazione

la funzione ha come parametri di default:

max_iter = 50, sono le partenze casuali massime
max_samples = 30, la dimensione del sottoinsieme di predittori utilizzati per creare le diverse interazioni possibili
improve = p-value = \(0.1\), pct = 0, la soglia per cui il t-test sul miglioramento è statisticamente significativo

set.seed(236)
ames_search <- fsa_two_way(ames_folds, ames_rec, lr_spec, multi_metric)

Risultati

## # A tibble: 31 x 4
##    var_1                var_2                                       pval   RMSE
##    <chr>                <chr>                                      <dbl>  <dbl>
##  1 Foundation: CBlock   Living Area                             0.00539  0.0513
##  2 Foundation: PConc    Living Area                             0.00828  0.0515
##  3 Building Type Duplex Living Area                             0.00178  0.0516
##  4 Living Area          MS SubClass: Duplex All Styles and Ages 0.00102  0.0516
##  5 Roof Style: Hip      Year Built                              0.0495   0.0518
##  6 Roof Style: Gable    Year Built                              0.000320 0.0518
##  7 Longitude (spline)   Neighborhood: Northridge Heights        0.00789  0.0518
##  8 Living Area          Neighborhood: Northridge Heights        0.0964   0.0519
##  9 Full Bath            Land Contour: HLS                       0.0767   0.0519
## 10 Foundation: CBlock   Lot Area                                0.0170   0.0519
## # ... with 21 more rows

L’output dell’algoritmo mostra le interazioni ordinate per il valore dell’RSME ottenuto dal modello con i main effect più l’interazione in questione. Inoltre è presente anche il valore del p-value ottenuto nel t-test per verificare se il miglioramento apportato dall’aggiunta dell’interazione al modello base si può ritenere significativo.

Per esempio sembra esserci una potenziale interazione tra il building type e la living area (p-value basso 0,00178).

Esistono altri modelli utili a scoprire le principali interazioni tra regressori basati sulla Multivariate adaptive regression splines (MARS) e sul modello Cubist.

In generale la Multivariate adaptive regression splines (MARS) è un modello di regressione non lineare per una variabile risposta di tipo continuo, essa ha come obiettivo la creazione di una hinge function utile a descrivere la relazione tra il predittore e la variabile dipendente. Tale funzione consente al modello di cercare relazioni non lineari. Nel contesto di ricerca dei principali termini di interazione, la MARS può essere utilizzata per cercare dei prodotti tra regressori utili a creare delle interazioni non lineari che isolano porzioni dello spazio delle variabili.

Cubist è un modello di regressione rule-based che costruisce un albero iniziale e lo scompone in un insieme di regole che vengono potate e qualora necessario eliminate. Per ciascuna regola viene definito un apposito modello lineare, creando un insieme di interazioni disgiunte. Le previsioni possono appartenere a più regole e, in questo caso, viene calcolata la media delle previsioni del modello lineare associato. Questa struttura è molto flessibile e ha una forte probabilità di individuare le interazioni locali all’interno dell’intero spazio dei predittori.

L’errore di previsone \(ErrF\) è definito come:

\[ErrF=E(MSE_{Te})=E \left[\frac{1}{n}\sum_{i=1}^{n} (Y^*_i-\hat f(x_i))^2\right] = \frac{1}{n}\sum_{i=1}^{n} E\left [(Y^*_i-\hat f(x_i))^2\right]\] E’ possibile distinguere tre tipologie di fonti che causano un errore di previsione:

errore irriducibile: dovuto alla presenza della componente stocastica contenuta nel modello.

Errore riducibile, scomponibile in:

Distorsione: esprime la quota di errore legata alla differenza media tra \(f\) e \(\hat f\)
Varianza: variabilità dello stimatore \(\hat f\)

Non esiste un modello che minimizza contemporaneamente la distorsione e la varianza, è necessario cercare un “compromesso”.

\[E[(Y^*_i-\hat Y^*_i)^2]=E[(f(x_i) + \epsilon^*_i - \hat f(x_i))^2]= \underbrace{E[(f(x_i) - \hat f(x_i))^2]}_{riducibile}+\underbrace{Var(\epsilon^*_i)}_{non\ riducibile} \] dove, \(Var(\epsilon^*_i)=\sigma^2\).

L’errore riducibile può essere scomposto nel \([Bias(\hat f(x_i))]^2\) e nella varianza \(Var( \hat f(x_i))\) dello stimatore \(\hat f\), rispettivamente:

\[E[ ( f(x_i) - \hat f(x_i) )^2 ] = \underbrace{[ E ( \hat f(x_i) - f(x_i) ) ]^2}_{[Bias( \hat f(x_i) ) ]^2 } + \underbrace{Var[\hat f(x_i)]}_{Variance[\hat f(x_i)]}\]

A partire da questa scomposizione l’errore di previsione \(ErrF\) si può scrivere come:

\[ErrF=\sigma^2+\underbrace{\frac{1}{n}\sum_{i=1}^{n} (E\hat f(x_i)-f(x_i))^2}_{Bias^2}+\underbrace{\frac{1}{n} \sum_{i=i}^{n} Var(\hat f(x_i) ) }_{Varianza}\]

Bias e varianza non possono essere minimizzate simultaneamente, qualora si fosse interessati a diminuire il bias si dovrà utilizzare un modello molto complesso ma in questo modo si otterrà un aumento della varianza. Esiste un trade-off tra bias e varianza:

Step di analisi (PSEUDO CODE)

Considerare il dataset ames come Population_data.
Dividere Population_data in Test_data e Training_data.
Costruire Population_model: addestrato su Population_data e applicato su Test_data.
Costruire il Mean_model, a partire da 20 campioni casuali estratti da Training_data è stato costruito un modello per ognuno di essi. Ciò permette di ottenere la previsione media dei modelli sul Test_data.
Calcolo Bias
Calcolo Varianza

library(caret)
library(glmnet)
library(tidymodels)
library(AmesHousing)
library(MASS)
library(ipred)
library(rpart)
library(randomForest)
library(reshape2)

Varianza

\[Variance (x) = E_{sample}[\hat f_{sample}(x) -\bar f(x)]^2\] Dove,

\(\bar f(x)\) modello medio, media delle previsioni ottenute a partire da 20 modelli costruiti su 20 campioni differenti.
\(\hat f_{sample}(x)\) valori previsti sul test set, stimati a partire da un modello costruito su un campione.

calculate_variance_of_model <- function(samplePredictions, y_test){
  predictions_mean_model <- colMeans(samplePredictions)
  colNames <- colnames(samplePredictions)
  variance = numeric()
  mse = numeric()
  i = 1
  for (colName in colNames) {
    variance[i] = mean(as.numeric(samplePredictions[,colName] - predictions_mean_model[i])^2)
    mse[i] = mean((samplePredictions[, colName] - as.numeric(y_test[i,]))^2)
    i=i+1
  }
  return(list(mean(variance), mean(mse)))
}

Bias

\[Bias (x) = (f(x) -\bar f(x))\]

Dove,

\(\bar f(x)\) modello medio, media delle previsioni ottenute a partire da 20 modelli costruiti su 20 campioni differenti.
\(f(x)\) modello generatore dei dati.

calculate_bias_of_model <- function(samplePredictions, y_hat_pop){
  predictions_mean_model <- colMeans(samplePredictions)
  return((mean(abs(predictions_mean_model-y_hat_pop)))^2)
}

Caricamento dataset

ames <- make_ames()

set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)

load("C:/Users/Stefano/Downloads/ames_h_stats (1).RData")

Creazione dei termini di interazione

int_vars <- 
  h_stats %>% 
  dplyr::filter(Estimator == "Bootstrap" & H > 0.001) %>% 
  pull(Predictor)


interactions <- t(combn(as.character(int_vars), 2))
colnames(interactions) <- c("var1", "var2")

interactions <- 
  interactions %>% 
  as_tibble() %>% 
  mutate(
    term = 
      paste0(
        "starts_with('",
        var1,
        "'):starts_with('",
        var2,
        "')"
      )
  ) %>% 
  pull(term) %>% 
  paste(collapse = "+")

interactions <- paste("~", interactions)
interactions <- as.formula(interactions)


int_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>%
  step_interact(interactions) %>% 
  step_zv(all_predictors()) %>%
  step_center(all_predictors()) %>%
  step_scale(all_predictors()) %>%
  step_bs(Longitude, Latitude, options = list(df = 5))

prep.int_rec <- prep(int_rec)
train_data <- bake(prep.int_rec, new_data = ames_train)
test_data <- bake(prep.int_rec, new_data = ames_test)
ames_data <- bake(prep.int_rec, new_data = ames)

#x_train <- train_data[,-8]
#y_train <- train_data[,8]

x_test <- test_data[,-8]
y_test <- test_data[,8]

#x_pop <- ames_data[,-8]
#y_pop <- ames_data[,8]

Calcolo modello popolazione

Considerato che la funzione generatrice dei dati è ignota, si assume che tale coincida con il modello stimato sulla totalità dei dati a disposizione.

In questa applicazione sono state considerate le seguenti funzioni generatrici:

TREE

fit.tree_population <- rpart(Sale_Price ~ ., ames_data)
phat.tree_population <- predict(fit.tree_population, newdata = x_test)

RANDOM FOREST

fit.rf <- randomForest(Sale_Price ~ ., data = ames_data, ntree = 50, importance = T)
phat.random_forest <- predict(fit.rf, newdata = x_test)

BAGGING

set.seed(123)
ames_bag1 <- bagging(
  formula = Sale_Price ~ .,
  data = ames_data,
  nbagg = 10,  
  coob = TRUE,
  control = rpart.control(minsplit = 2, cp = 0)
)
y_hat_pop_bag<- predict(ames_bag1, newdata = x_test)

Stima dei modelli

Al fine di stimare i modelli sui diversi campioni, si sono implementate le seguenti funzioni che permettono di addestrare un modello e ottenere le previsioni sul test set.

TREE

samplePredForDecisionTree <- function(campione){
  sample_Tree_Model <- rpart(Sale_Price ~ ., campione)
  return(predict(sample_Tree_Model, newdata = x_test))
}

RANDOM FOREST

samplePredForRandomForest <- function(campione){
  sample_RF_Model <- randomForest(Sale_Price ~ ., data = campione, ntree = 50, importance = T)
  return(predict(sample_RF_Model, newdata = x_test))
}

BAGGING

samplePredBagging <- function(campione){
  set.seed(123)
  library(caret)
  library(ipred)
  # train bagged model
  ames_bag <- bagging(
    formula = Sale_Price ~ .,
    data = campione,
    nbagg = 10,  
    coob = TRUE,
    control = rpart.control(minsplit = 2, cp = 0)
  )
  return(predict(ames_bag, newdata = x_test))
}

Calcolo bias & varianza

In questa applicazione vengono presi in considerazione 20 modelli per ciascuna ampiezza campionaria. Le ampiezze campionarie di riferimento sono da 100, 300, 500, 700, 900, 1000, 1200 osservazioni.

Con riferimento alla prima ampiezza campionaria formata da 100 osservazioni, si procede costruendo 20 campioni, ciascuno dei quali è ottenuto tramite un campionamento casuale senza reinserimento. Per ogni campione viene addestrato un generico modello, le cui previsioni sono stimate sul test set. Ciò permette di calcolare empiricamente il bias e la varianza del modello.

Si ripete questa procedura per le diverse ampiezze campionarie, in questo modo è possibile valutare l’andamento del bias e varianza all’aumentare della dimensione campionaria.

nModel <-20
noOfSamples = c(100,300,500,700,900,1000,1200)
samplePredictionsTree <- data.frame()
samplePredictionsRF <- data.frame()
samplePredictionsBag <- data.frame()
risultato_finale = list()
risultato_finale_rf = list()
risultato_finale_bg = list()

for (numerosita in noOfSamples) {
  for (i in 1:nModel) {
    campione = train_data[sample(nrow(ames_train),numerosita, replace = FALSE),]
    samplePredictionsTree_temp = samplePredForDecisionTree(campione)
    samplePredictionsTree <- rbind(samplePredictionsTree,samplePredictionsTree_temp)
    samplePredictionsRF_temp = samplePredForRandomForest(campione)
    samplePredictionsRF <- rbind(samplePredictionsRF,samplePredictionsRF_temp)
    samplePredictionsBag_temp = samplePredBagging(campione)
    samplePredictionsBag <- rbind(samplePredictionsBag,samplePredictionsBag_temp)
  }
  var_model = calculate_variance_of_model(samplePredictionsTree, y_test)
  bias_model = calculate_bias_of_model(samplePredictionsTree, phat.tree_population)
  var_model_rf = calculate_variance_of_model(samplePredictionsRF, y_test)
  bias_model_rf = calculate_bias_of_model(samplePredictionsRF, phat.random_forest)
  var_model_bg = calculate_variance_of_model(samplePredictionsBag, y_test)
  bias_model_bg = calculate_bias_of_model(samplePredictionsBag, y_hat_pop_bag)
  
  risultato_finale_bg[[numerosita]] <- c(numerosita,var_model_bg[[1]], var_model_bg[[2]],bias_model_bg)
  risultato_finale_rf[[numerosita]] <- c(numerosita,var_model_rf[[1]], var_model_rf[[2]],bias_model_rf)
  risultato_finale[[numerosita]] <- c(numerosita,var_model[[1]], var_model[[2]],bias_model)
}

Le seguenti rappresentazioni grafiche mostrano l’andamento di \(bias^2\), \(varianza\) e \(MSE\) al variare dell’ampiezza campionaria. All’aumentare dell’ampiezza campionaria aumenta la \(varianza\), ma si riduce il \(bias^2\), non esiste un modello che minimizza contemporaneamente la distorsione e la varianza, è necessario determinare un “compromesso”.

Decision Tree

risultato_finale <- do.call("rbind",risultato_finale)
colnames(risultato_finale) <- c("numerosita","varianza","mse", "bias")
res_tree <- as_tibble(risultato_finale)
res_tree$method <-"Decision Tree"
res_tree

## # A tibble: 7 x 5
##   numerosita varianza    mse    bias method       
##        <dbl>    <dbl>  <dbl>   <dbl> <chr>        
## 1        100 0.000705 0.0144 0.00536 Decision Tree
## 2        300 0.00272  0.0136 0.00280 Decision Tree
## 3        500 0.00345  0.0128 0.00217 Decision Tree
## 4        700 0.00344  0.0121 0.00157 Decision Tree
## 5        900 0.00317  0.0118 0.00134 Decision Tree
## 6       1000 0.00308  0.0115 0.00135 Decision Tree
## 7       1200 0.00292  0.0112 0.00131 Decision Tree

require(gridExtra)
plot1 <- ggplot(res_tree, aes(x=numerosita, y=bias))+
  geom_line(stat = "identity", col="blue")+
  ggtitle("Bias")

plot2 <- ggplot(res_tree, aes(x=numerosita, y=varianza))+
  geom_line(stat = "identity", col="red")+
  ggtitle("Varianza")

plot3 <-ggplot(res_tree, aes(x=numerosita, y=mse))+
  geom_line(stat = "identity", col="green")+
  ggtitle("MSE")

grid.arrange(plot1, plot2, plot3, ncol=2)

Random Forest

risultato_finale_rf <- do.call("rbind",risultato_finale_rf)
colnames(risultato_finale_rf) <- c("numerosita","varianza","mse", "bias")
res_rf <- as_tibble(risultato_finale_rf)
res_rf$method <-"Random Forest"
res_rf

## # A tibble: 7 x 5
##   numerosita varianza     mse    bias method       
##        <dbl>    <dbl>   <dbl>   <dbl> <chr>        
## 1        100 0.000195 0.00830 0.00257 Random Forest
## 2        300 0.000511 0.00819 0.00223 Random Forest
## 3        500 0.000853 0.00762 0.00184 Random Forest
## 4        700 0.000913 0.00710 0.00157 Random Forest
## 5        900 0.000897 0.00681 0.00144 Random Forest
## 6       1000 0.000913 0.00673 0.00137 Random Forest
## 7       1200 0.000878 0.00654 0.00130 Random Forest

plot1 <- ggplot(res_rf, aes(x=numerosita, y=bias))+
  geom_line(stat = "identity", col="blue")+
  ggtitle("Bias")

plot2 <- ggplot(res_rf, aes(x=numerosita, y=varianza))+
  geom_line(stat = "identity", col="red")+
  ggtitle("Varianza")

plot3 <- ggplot(res_rf, aes(x=numerosita, y=mse))+
  geom_line(stat = "identity", col="green")+
  ggtitle("MSE")

grid.arrange(plot1, plot2, plot3, ncol=2)

Bagging

risultato_finale_bg <- do.call("rbind",risultato_finale_bg) 
colnames(risultato_finale_bg) <- c("numerosita","varianza","mse", "bias")
res_bg <- as_tibble(risultato_finale_bg)
res_bg$method <-"Bagging"
res_bg

## # A tibble: 7 x 5
##   numerosita varianza     mse    bias method 
##        <dbl>    <dbl>   <dbl>   <dbl> <chr>  
## 1        100 0.000356 0.00975 0.00354 Bagging
## 2        300 0.00135  0.0101  0.00301 Bagging
## 3        500 0.00195  0.00944 0.00250 Bagging
## 4        700 0.00201  0.00875 0.00213 Bagging
## 5        900 0.00194  0.00845 0.00201 Bagging
## 6       1000 0.00198  0.00833 0.00191 Bagging
## 7       1200 0.00189  0.00803 0.00182 Bagging

plot1 <- ggplot(res_bg, aes(x=numerosita, y=bias))+
  geom_line(stat = "identity", col="blue")+
  ggtitle("Bias")

plot2 <- ggplot(res_bg, aes(x=numerosita, y=varianza))+
  geom_line(stat = "identity", col="red")+
  ggtitle("Varianza")

plot3 <- ggplot(res_bg, aes(x=numerosita, y=mse))+
  geom_line(stat = "identity", col="green")+
  ggtitle("MSE")

grid.arrange(plot1, plot2, plot3, ncol=2)

La seguente rappresentazione confronta \(MSE\), \(bias^2\) e \(varianza\) per i tre modelli considerati fissata un’ampiezza campionaria pari a 1200 osservazioni. Dai risultati emerge che il random forest ha varianza minore, infatti a differenza del bagging, per ciascun nodo anziché esplorare tutti i possibili predittori, viene preso in considerazione solamente un sottoinsieme costituito da \(m<p\) predittori.

res_tot <- rbind(res_tree[7,], res_rf[7,],res_bg[7,])

ggplot(melt(res_tot[,2:5], id.vars='method'),aes(method,value,fill=variable))+
     geom_bar(stat="identity",position="dodge")

In generale, in termini di errore quadratico medio il random forest risulta essere il modello migliore.

NB: Il Bagging è un caso particolare di Random Forest, infatti ponendo \(m = p\) si ottengono gli stessi risultati. Nel Random Forest le righe vengono ricampionate nello stesso modo del bagging (ricampionamento con reinserimento) ma viene considerato solamente un sottoinsieme dei predittori. Il difetto principale del Bagging deriva dal limite insito del bootstrap: ovvero, i training set possono essere molto simili tra loro portando così ad ottenere degli alberi molto correlati. Il Random Forest si pone l’obiettivo di rendere gli alberi costruiti dai campioni bootstrap più diversi possibili.

Analisi per modello glmnet e lasso

Carico dati e splitting in training e test set

library(caret)
library(glmnet)
library(tidymodels)
library(AmesHousing)
library(hdi)
load("C:/Users/Stefano/Downloads/ames_h_stats (1).RData")

ames <- make_ames()

set.seed(955)
ames_split <- initial_split(ames)
ames_train <- training(ames_split)
ames_test <- testing(ames_split)

set.seed(24873)
ames_folds <- vfold_cv(ames_train)
ames_ind <- rsample2caret(ames_folds)

Costruzione del modello senza interazioni:

Si costruisce un modello con variabile target Sale_Price e con 18 variabili esplicative: Bldg_Type, Neighborhood, Year_Built, Gr_Liv_Area, Full_Bath, Year_Sold, Lot_Area, Central_Air, Longitude, Latitude, MS_SubClass, Alley, Lot_Frontage, Pool_Area, Garage_Finish, Foundation, Land_Contour, Roof_Style. La variabile target viene sottoposta a trasformazione logaritmica, mentre le variabili Lot_Area, Gr_Liv_Area, Lot_Frontage vengono normalizzate con una trasformazione di BoxCox. La dimensionalità della variabile Neighborhood viene ridotta includendo in un’unica classe “other” tutte le osservazioni con frequenza relativa <0.05. Tutte le variabili nominali vengono trasformate in variabili dummy e tutte le variabili esplicative vengono normalizzate o eliminate se hanno un solo valore identico in ogni osservazione.

main_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>%
  step_zv(all_predictors()) %>%
  step_center(all_predictors()) %>%
  step_scale(all_predictors()) %>%
  step_bs(Longitude, Latitude, options = list(df = 5))
main_rec

## Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         18
## 
## Operations:
## 
## Log transformation on Sale_Price
## Box-Cox transformation on Lot_Area, Gr_Liv_Area, Lot_Frontage
## Collapsing factor levels for Neighborhood
## Dummy variables from all_nominal()
## Zero variance filter on all_predictors()
## Centering for all_predictors()
## Scaling for all_predictors()
## B-Splines on Longitude, Latitude

Selezione e costruzione dei termini di interazione rilevanti in base alla statistica H:

Considerata la distribuzione della statistica H ottenuta a partire dal metodo bootstrap, si sono identificati i potenziali predittori fissando un cutoff pari a 0.001.

int_vars <- 
  h_stats %>% 
  dplyr::filter(Estimator == "Bootstrap" & H > 0.001) %>% 
  pull(Predictor)

Si procede creando un dataframe contenente tutte le possibile coppie di interazioni \(\frac{p(p−1)}{2}\) tra i regressori aventi un valore della statistica H maggiore della soglia.

interactions <- t(combn(as.character(int_vars), 2))
colnames(interactions) <- c("var1", "var2")
interactions

##       var1            var2          
##  [1,] "Garage_Finish" "Roof_Style"  
##  [2,] "Garage_Finish" "Central_Air" 
##  [3,] "Garage_Finish" "Year_Built"  
##  [4,] "Garage_Finish" "Neighborhood"
##  [5,] "Garage_Finish" "Gr_Liv_Area" 
##  [6,] "Roof_Style"    "Central_Air" 
##  [7,] "Roof_Style"    "Year_Built"  
##  [8,] "Roof_Style"    "Neighborhood"
##  [9,] "Roof_Style"    "Gr_Liv_Area" 
## [10,] "Central_Air"   "Year_Built"  
## [11,] "Central_Air"   "Neighborhood"
## [12,] "Central_Air"   "Gr_Liv_Area" 
## [13,] "Year_Built"    "Neighborhood"
## [14,] "Year_Built"    "Gr_Liv_Area" 
## [15,] "Neighborhood"  "Gr_Liv_Area"

interactions <- 
  interactions %>% 
  as_tibble() %>% 
  mutate(
    term = 
      paste0(
        "starts_with('",
        var1,
        "'):starts_with('",
        var2,
        "')"
      )
  ) %>% 
  pull(term) %>% 
  paste(collapse = "+")

interactions <- paste("~", interactions)
interactions <- as.formula(interactions)

Costruzione del modello glmnet con interazioni:

int_rec <-
  recipe(Sale_Price ~ Bldg_Type + Neighborhood + Year_Built +
           Gr_Liv_Area + Full_Bath + Year_Sold + Lot_Area +
           Central_Air + Longitude + Latitude + MS_SubClass +
           Alley + Lot_Frontage + Pool_Area + Garage_Finish + 
           Foundation + Land_Contour + Roof_Style,
         data = ames_train) %>%
  step_log(Sale_Price, base = 10) %>%
  step_BoxCox(Lot_Area, Gr_Liv_Area, Lot_Frontage) %>%
  step_other(Neighborhood, threshold = 0.05) %>% 
  step_dummy(all_nominal()) %>%
  step_interact(interactions) %>% 
  step_zv(all_predictors()) %>%
  step_center(all_predictors()) %>%
  step_scale(all_predictors()) %>%
  step_bs(Longitude, Latitude, options = list(df = 5))

Parameter tuning del training:

La ricerca dei valori ottimi di \(\alpha\) e \(\lambda\) avviene in una griglia contenente i possibili valori per tali parametri.

ctrl <- 
  trainControl(
    method = "cv",
    index = ames_ind$index,
    indexOut = ames_ind$indexOut
  )

glmn_grid <- expand.grid(alpha = seq(.2, 1, by = .2), lambda = 10^seq(-4, -1, by = 0.1))

Training del modello glmnet con interazioni:

main_glmn_h <- 
  train(main_rec,
        data = ames_train, 
        method = "glmnet",
        tuneGrid = glmn_grid,
        trControl = ctrl
  )

Training del modello glmnet senza interazioni:

int_glmn_h <- 
  train(int_rec,
        data = ames_train, 
        method = "glmnet",
        tuneGrid = glmn_grid,
        trControl = ctrl
        )

int_glmn_h$bestTune

##     alpha      lambda
## 136     1 0.001258925

Plot parametri di tuning di int_glmn_h:

p <- 
  ggplot(int_glmn_h) + 
  scale_x_log10() + 
  theme_bw() + 
  theme(legend.position = "top")
p

Questo grafico mostra l’andamento dell’RMSE al variare dei parametri \(\lambda\) e \(\alpha\). Ogni curva corrisponde ad un diverso valore di \(\alpha\) e rappresenta le performance del modello stimato in funzione del parametro di penalizzazione \(\lambda\). In questo modo è possibile individuare i valori dei parametri di tuning che minimizzano l’RMSE.

Scelta di lambda con cross-validation

Il nostro obiettivo è quello di selezionare tra i 157 regressori quelli che risultano maggiormente implicati nello spiegare la variazione della variabile dipendente. La ricerca di lambda ottima avviene nella griglia di valori specificata all’interno del comando cv.glmnet.

lasso = cv.glmnet(as.matrix(train_data[,-8]), 
                  train_data$Sale_Price, alpha=1, 
                  lambda = c(10^seq(-4, -1, by = 0.1)))
plot(lasso, ylim = c(0,0.04))

In questo grafico si evince il lambda path: andamento della curva dell’errore quadratico medio di previsione dato dalla cross validation. Per ciascun valore di lambda la media dell’errore di cross validation è rappresentata dal punto rosso le barre rappresentano gli standard error della media in corrispodenza di quel valore di lambda.

## 
## Call:  cv.glmnet(x = as.matrix(train_data[, -8]), y = train_data$Sale_Price,      lambda = c(10^seq(-4, -1, by = 0.1)), alpha = 1) 
## 
## Measure: Mean-Squared Error 
## 
##       Lambda Index  Measure        SE Nonzero
## min 0.001259    20 0.005689 0.0005258      67
## 1se 0.005012    14 0.006049 0.0005793      32

Il valore di \(\lambda\) che minimizza la curva di cross validation permette di stabilire che il numero di predittori con coefficente di regressione diverso da 0 è 67.

NB.
É possibile “tollerare” un piccolo incremento nell errore guadagnandoci dal punto di vista interpretativo; in questo modo scelgo un nuovo valore di \(\lambda_{1se}\) che permette di ottenere un errore medio di previsione maggiore ma al contempo restituisce un modello più interpretabile.

Performance dei modelli sul test set

Per valutare l’effetto delle interazioni sulla capacità di prevedere la varibile dipendente, si confronta l’errore quadratico medio dei due modelli sul test set.

yhats_main = predict(main_glmn_h, ames_test)
yhats_int = predict(int_glmn_h, ames_test)

## [1] "Modello main:  0.0053"

## [1] "Modello main + interazioni:  0.0054"

In termini di errore quadratico medio i due modelli risultano pressoché uguali, nonostante questo si predilige il modello senza interazioni perchè meno complesso e più facilmente interpretabile.

Trade-off bias & varianza

Calcolo modello popolazione

Considerando che la funzione generatrice dei dati è ignota, si assume che questa coincida con il modello stimato sulla totalità dei dati a nostra disposizione.

In questa applicazione è stata considerata la seguente funzione generatrice, che coincide con il modello lasso caratterizzato dal valore di \(\lambda\) che minimizza l’errore quadratico medio.

pop_lasso_model <- glmnet(ames_data[,-8], ames_data$Sale_Price, alpha = 1, lambda = 0.001258925)
y_hat_pop_lasso <- predict(pop_lasso_model, newx = as.matrix(x_test))

Stima del modello

La seguente funzione consente di addestrare il modello e ottenere le previsioni sul test set.

samplePredForLasso <- function(campione, lambda){
  sample_Lasso_Model <- glmnet(campione[,-8], campione$Sale_Price, alpha = 1, lambda = lambda)
  return(predict(sample_Lasso_Model, newx = as.matrix(x_test)))
}

Calcolo di bias & varianza

In questa applicazione vengono stimati 30 modelli con un’ampiezza campionaria pari a 2000, ciascuno dei quali è ottenuto tramite un campionamento casuale con reinserimento.

Fissato un valore di \(\lambda\), per ogni campione viene stimato un modello lasso e si effettuano previsioni sul test set. Ciò permette di calcolare empiricamente il bias e la varianza del modello.

Si ripete questa procedura per i diversi valori di \(\lambda\), in questo modo è possibile valutare l’andamento del bias e varianza all’aumentare della dimensione campionaria.

nModel <-30
samplePredictionsLasso <- data.frame()
risultato_finale = list()
numerosita=2000
lambdas = c(0.001258925, 0.01, 0.05,0.1, 0.15, 0.2, 0.25, 0.30, 0.35, 0.40)
j=1
for (lambda in lambdas) {
  for (i in 1:nModel) {
    campione = train_data[sample(nrow(ames_train),numerosita, replace = TRUE),]
    samplePredictionsLasso_temp = samplePredForLasso(campione, lambda)
    samplePredictionsLasso <- rbind(samplePredictionsLasso,samplePredictionsLasso_temp)
  }
  var_model = calculate_varaince_of_model(samplePredictionsLasso, y_test)#aggiungere mse
  bias_model = calculate_bias_of_model(samplePredictionsLasso, y_hat_pop_lasso)#aggiungere mse
  
 
  risultato_finale[[j]] <- c(lambda,var_model[[1]], var_model[[2]],bias_model)
  j=j+1 
}

risultato_finale <- do.call("rbind",risultato_finale)
colnames(risultato_finale) <- c("lambda","varianza","mse", "bias")
risultato_finale <- as_tibble(risultato_finale)
risultato_finale

## # A tibble: 10 x 4
##     lambda varianza    mse   bias
##      <dbl>    <dbl>  <dbl>  <dbl>
##  1 0.00126  0.0254  0.0386 0.0168
##  2 0.01     0.0229  0.0363 0.0168
##  3 0.05     0.0183  0.0316 0.0168
##  4 0.1      0.0142  0.0273 0.0168
##  5 0.15     0.0114  0.0244 0.0168
##  6 0.2      0.00949 0.0224 0.0168
##  7 0.25     0.00814 0.0209 0.0168
##  8 0.3      0.00712 0.0198 0.0168
##  9 0.35     0.00634 0.0190 0.0169
## 10 0.4      0.00570 0.0183 0.0169

ggplot(risultato_finale, aes(x=lambda, y=bias))+
  geom_line(stat = "identity", col="blue")+
  ggtitle("Bias")

ggplot(risultato_finale, aes(x=lambda, y=varianza))+
  geom_line(stat = "identity", col="red")+
  ggtitle("Varianza")

ggplot(risultato_finale, aes(x=lambda, y=mse))+
  geom_line(stat = "identity", col="green")+
  ggtitle("MSE")

Da questi grafici si evince la presenza di un trade-off tra \(distorsione^2\) e varianza, in quanto all’aumentare della penalizzazione diminuisce la varianza, ma aumenta il \(bias^2\); questo perchè il modello lasso esclude un numero sempre maggiore di predittori.

Detecting interaction effect

Cattaneo Andrea, Cosentino Alessandro, Marzorati Stefano

9/12/2021

ESEMPIO DI INTERAZIONE TRA UNA VARIABILE QUANTITATIVA E UNA VARIABILE QUALITATIVA

TIPOLOGIE DI INTERAZIONE

PRINCIPI PER LA RICERCA DELLE INTERAZIONI

Analisi esplorativa

Qual è il costo di una casa?

Quando sono state costruite le case?

Quando sono state vendute le case?

In che quartiere sono situate le case?

Quanto sono grandi le case?

Analisi Univariata

Correlation matrix

The Brute-Force Approach to Identifying Predictive Interactions

Simple Screening

Esempio dataset Ames

Penalized Regression

Glmnet regression

Esempio con dati ames

Recupero dati e creazione modello main

Creazione interazioni e modelli main + interazioni

Tuning dei modelli

Modello main

Modello main + interazioni

Risultato main + interazioni

Creazione modello two-way

Tree-Based Methods

Esempio dati ames

The Feasible Solution Algorithm

Esempio con dati ames

Recupero dati e creazione modello base

Inizializzazione parametri

Esecuzione dell’algoritmo fsa

Risultati

Sviluppi futuri

Analisi trade-off Bias & Varianza dei modelli con termini di interazione

Step di analisi (PSEUDO CODE)

Varianza

Bias

Caricamento dataset

Creazione dei termini di interazione

Calcolo modello popolazione

Stima dei modelli

Calcolo bias & varianza

Decision Tree

Random Forest

Bagging

Analisi per modello glmnet e lasso

Carico dati e splitting in training e test set

Costruzione del modello senza interazioni:

Selezione e costruzione dei termini di interazione rilevanti in base alla statistica H:

Costruzione del modello glmnet con interazioni:

Parameter tuning del training:

Training del modello glmnet con interazioni:

Training del modello glmnet senza interazioni:

Plot parametri di tuning di int_glmn_h:

Scelta di lambda con cross-validation

Performance dei modelli sul test set

Trade-off bias & varianza

Calcolo modello popolazione

Stima del modello

Calcolo di bias & varianza